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Abstract. In this paper we study a model of geometry of vision due to Petitot, Citti and Sarti. 
One of the main features of this model is that the primary visual cortex VI lifts an image from ]R2 to 
the bundle of directions of the plane. Neurons are grouped into orientation columns, each of them 
corresponding to a point of this bundle. 

In this model a corrupted image is reconstructed by minimizing the energy necessary for the 
activation of the orientation columns corresponding to regions in which the image is corrupted. The 
minimization process intrinsically defines an hypoelliptic heat equation on the bundle of directions 
of the plane. 

In the original model, directions are considered both with and without orientation, giving rise 
respectively to a problem on the group of rototranslations of the plane SE{2) or on the projective 
tangent bundle of the plane PTR^ . 

We provide a mathematical proof of several important facts for this model. We first prove that 
the model is mathematically consistent only if directions are considered without orientation. We then 
prove that the convolution of a L^(R^, R) function (e.g. an image) with a 2-D Gaussian is generically 
a Morse function. This fact is important since the lift of Morse functions to PTR^ is defined on a 
smooth manifold. We then provide the explicit expression of the hypoelliptic heat kernel on PTR-^ 
in terms of Mathieu functions. 

Finally, we present the main ideas of an algorithm which allows to perform image reconstruction 
on real non-academic images. A very interesting point is that this algorithm is massively parallelizable 
and needs no information on where the image is corrupted. 

Keywords: sub-Riemannian geometry, image reconstruction, hypoelliptic diffu- 
sion 

1. Introduction. In this paper we study a model of geometry of vision due to 
Petitot, Citti and Sarti. The main reference for the model is the paper [15]. Its first 
version can be found in [37, 39]. This model was also studied by the authors of the 
present paper in [10[, by Hladky and Pauls [24[ and, independently, by Duits et al. 
in a series of papers mostly for contour completion [17] and contour enhancement 
[18, 19[. This model has been called the pinwheel model by Petitot himself, see [40]. 
See also ]38, 45] and references therein. 

To start with, assume that a grey-level image is represented by a function I G 
L^{'D, R), where V is an open bounded domain of R^. The algorithm that we present 
here is based on three crucial ideas coming from neurophysiology: 

1. It is widely accepted that the retina approximately smoothes the images by 
making the convolution with a Gaussian function (see for instance ]28, 33, 36] 
and references therein), equivalently solving a certain isotropic heat equation. 
Moreover, smoothing by the same technique is a widely used method in image 
processing. Then, it is an interesting question in itself to understand generic 
properties of these smoothed images. Our first result (proved in Appendix 
A) is that, given G{(Tx,<Ty) the two dimensional Gaussian centered in (0,0) 
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with standard deviations a^jO-y > 0, then the smoothed image 



is generically a Morse function (i.e. a smooth function having as critical 
points only non-degenerate maxima, minima and saddles). This has interest- 
ing consequences, as explained in the following. 

Remark 1. In several applications, the convolution is made with a Gaus- 
sian of small standard deviations. Equivalently, the smoothed image can be 
obtained as the solution of an isotropic heat equation with small final time. 
Remark 2. These results can be generalized to non-Gaussian filters and 
even to non-linear smoothing processes. See for instance [16] for some of 
these generalizations. 
2. The primary visual cortex VI lifts the image from to the bundle 
of directions of the plane PTM.^. 

In a simplified model"'^ (see [15] and [38, p. 79]), neurons of VI are grouped into 
orientation columns, each of them being sensitive to visual stimuli at a given 
point a of the retina and for a given direction p on it. The retina is modeled by 
the real plane, i.e. each point is represented by a € M^, while the directions 
at a given point arc modeled by the projective line, i.e. p € P^. Hence, 
the primary visual cortex VI is modeled by the so called projective tangent 
bundle PTR^ := x P^. From a neurological point of view, orientation 
columns are in turn grouped into hypercolumns, each of them being sensitive 
to stimuli at a given point a with any direction. In the same hypercolumn, 
relative to a point a of the plane, we also find neurons that are sensitive to 
other stimuli properties, like colors. In this paper, we focus only on directions 
and therefore each hypercolumn is represented by a fiber P^ of the bundle 
PTR2. See Figure 1.1. 

This space has the topology of x P^ (it is a trivial bundle) and its points 
are triples {x, y, 0), where {x, y) gM.'^, 9 G K/(7rZ). 

The smoothed image / : — >^ R is lifted to a a function / defined as follows: 



It follows that / has support on a set 5/ C PTM^. The following fact consti- 
tutes our second result. If / is a Morse function (which happens generically 
due to the smoothing of the retina as explained above), then Sf is an embed- 
ded surface in PTR^, see Proposition 19. 
3. If the image is corrupted or missing on a set C I? (i.e. if I is defined on 
V then the reconstruction in is made by minimizing a given cost. 

This cost represents the energy that the primary visual cortex should spend 
in order to excite orientation columns which corresponds to points in fl and 
hence that are not directly excited by the image. An orientation column is 
easily excited if it is close to another (already activated) orientation column 
sensitive to a similar direction in a close position (i.e. if the two are close in 



^For example, in this model we do not tcike into account the fact that the continuous space of 
stimuli is implemented via a discrete set of neurons. 




f{x, y) if 9 is the direction of the level set of /, 



otherwise. 



PTM2). 
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Figure 1.1. A scheme of the primary visual cortex VI. 



When the image to be reconstructed is a curve, this gives rise to a sub- 
Riemannian problem (i.e. an optimal control problem linear in the control 
and with quadratic cost) on PTM^, which we briefly discuss in Sections 2.3, 
2.4, 2.5. 2 

When the image is not just a curve, the reconstruction is made by considering 
the difFusion process naturally associated with the sub-Riemannian problem 
on PTM^ (described by an hypoelliptic heat equation). Such a reconstruc- 
tion makes use of the function / as initial condition in a suitable way. The 
reconstructed image is then obtained by projecting the result of the diffusion 
from PTM? to R^. 

In this paper we study this model, providing a mathematical proof of several key 
facts and adding certain important details with respect to its original version given 
in [15, 45]. The main improvements are the following: 

A) As already mentioned, we start with any function I G Lp'i'D, M) and we prove that 
after convolution with a Gaussian of standard deviations'^ Ux = <^y, generically, 
we are left with a Morse function / e L2(]R2^r) n C°°(R2,M) (see Appendix 
A). This smoothing process is important to guarantee certain regularity of the 
domain of definition of the lifted function /. 

B) Our definition of the lift is suitable to all smooth functions, since we don't require 
conditions like nondegenerate gradient (as in [15]) or more complicated condition 
on the so called non-Legendrian solitary points (as in [24, Thm 1.6]). 

C) Recall that PTR^ can be seen as the quotient of the group of rototranslations 
of the plane SE{2) ~ R'^ x S'^ by Z2, where the quotient is the identification of 



■^In particular this sub-Riemannian structure has an underlying contact structure. To our knowl- 
edge, the first time in which the visual cortex was modeled as a contact structure was in [25]. 
•^We fix cTx = cTy to guarantee invariance by rototranslations of the algorithm. 
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antipodal points in S^. In the first version of this model [15] the image is lifted on 
SE{2) (i.e. directions arc considered with orientation), while in the seeond one 
[45] it is lifted on PTR^ (i.e. directions are considered without orientation). The 
next contribution of our paper is to show that the problem of reconstruction of 
images for smooth functions is well posed on PTR^ while it is not on SE{2). First, 
on PTR^ the lift is unique, while on SE{2) it is not, since level sets of the image 
are not oriented curves. Second, the problem on SE{2) cannot be interpreted as 
a problem of reconstruction of contours (see Remarks 4, 13 and [10]). Third, as 
proved in Proposition 19, the domain 5*/ of the lift of a Morse function / is much 
more natural on PTR^ than on SE{2). On PTR^ it is a manifold, while on SE{2) 
it is a manifold with a boundary (for a continuous choice of the orientation of the 
level sets of /). The boundary appears on minima, maxima and saddles of /. In 
the diffusion process, starting with an initial condition which is concentrated on 
a manifold is much more natural than starting with an initial condition which is 
concentrated on a manifold with a boundary. 

D) We show that the sub-Riemannian structure over PTM? is not trivializable, which 
means that it cannot be specified by a single global orthonormal frame as in [45]. 
For a detailed discussion of this issue see Remark 6 and [10]. 

E) We give the expression of the hypoelliptic heat kernel over PTR^, while, pre- 
viously, it was known only on SE{2) (see [4] and [17, 18], where it was found 
independently). 

F) We provide an effective algorithm for image reconstruction that looks unexpect- 
edly efficient on real non-academic examples as shown in Section 3.3. Moreover, 
our algorithm has the good feature to be massively parallelizable (see Section 
3). This is just the materialization of the classical fact that the noncommutative 
Fourier transform disintegrates the regular representation over SE{2). More- 
over the algorithm does not need the information of where the original image is 
corrupted. 

Other numerical methods to compute hypoelliptic diffusion on SE{2) for image 
processing have been developed. For instance: group convolution methods (see [14, 
18, 20]) finite differences [15, 45, 21]^ and finite elements methods. (See [18, 19]. 
These last works are related to the noncommutative Fourier transform on SE{2) and 
are extensions of the works by August [7].) Most of these works are about contour 
enhancement. 

Remark 3. Notice that from the very beginning of the algorithm, we deal with 
the intensity of the image. Other related algorithms [15, Sec. 3.3], [24] are instead 
composed of two reconstruction steps. After the lift of the image, these algorithms 
have to deal with a surface in SE{2) or PTR^ with a hole, corresponding to the 
corrupted part. The first reconstruction step is thus to fill the hole as a surface, 
without considering the intensity of the image. The second reconstruction step is 
then to put the intensity on the reconstructed part. See Remarks 15 and 22. 

The results of our algorithm can be compared to the ones coming from psycho- 
logical experiments. Moreover, they can be useful to reconstruct the geometry of an 
image, as a preliminary step of exemplar-based methods (see [12]). 

Notice that an alternative technique of image processing (in particular for contour 
completion) based on physiological models of the visual cortex has been proposed by 



^Notice that classical finite difference methods "hardly works" to compute hypoelliptic diffusion. 
This is due to the diffusion at different scales on different directions as a consequence of the non- 
ellipticity of the diffusion operator. 
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Mumford [32], then developed in [17, 18, 19]. In these models, contours (or level sets of 
an image) are considered with orientation and a non-isotropic diffusion is associated to 
an optimal control problem with drift having elastica curves as solutions. We briefly 
compare the method presented in this paper with the one by Mumford in Section 
2.8.1. 

The structure of the paper is the following. In Section 2 we present in detail 
the sub-Riemannian structure defined on PTM^. We then define the corresponding 
hypocUiptic diffusion which is one of the main tools used in the algorithm of image 
reconstruction and we find explicitly the corresponding kernel on PTM."^. At the end, 
we present in detail the mathematical algorithm. 

Section 3 is devoted to the discussion about the numerical integration of the 
hypoelliptic evolution and to the presentation of samples. 

Appendix A is devoted to the detailed proof that, generically, the convolution of 
a function over a bounded domain 2? C with a Gaussian G is a Morse function. 
In particular, we prove that the set of functions I £ L'^iP) whose convolution with a 
Gaussian is a Morse function in L'^{T>) is residual (i.e. it is a countable intersection 
of open and dense sets) . We then prove that the set of functions I <E (V) such that 
X*G restricted to a compact if C is a Morse function is open and dense. Notice 
again that the proof can be adapted to any reasonable smoothing process, not only 
Gaussian. 

2. The mathematical model and the algorithm. 

2.1. Reconstruction of a curve. In this section we briefly describe an algo- 
rithm to reconstruct interrupted planar curves. The main interest of this section is 
the definition of the sub-Riemannian structure over PTM?, from which we are going 
to define the sub-elliptic diffusion equation. The main reference for this algorithm is 
[15[, where the lift of a planar curve was defined on SE{2) rather than on PT^. 

Consider a smooth function 70 : [0,6] U [c, (i] — >■ (with a < h < c < d) 
representing a curve that is partially hidden or deleted in (6, c). We assume that 
starting and ending points never coincide, i.e. 70 (^) 7^ 7o(c), and that initial and final 
velocities 7(6) and 7(c) are well defined and nonvanishing. 

We want to find a curve 7 : [6, c] — > that completes 70 in the deleted part and 
that minimizes a cost depending both on the length and on the curvature of 7. 
Recall that 

^ xy - iji 
^ (x2 + y2)3/2 

where {x, y) are the components of 7. 

The fact that 7 completes 70 means that 7(6) = 7o(^), 7(c) = 70(c). It is also 
reasonable to require that the directions of tangent vectors coincide, i.e. 7(6) w 
io{h), 7(c) « 70(c) where 

vi « V2 if it exists a G M \ {0} such that Vi = av2- (2-1) 

Remark 4. Notice that we have required boundary conditions on initial and final 
directions without orientation. The problem above can also be formulated requiring 
boundary conditions with orientation, i.e. substituting in (2.1) the condition a € M"*". 
However, this choice docs not guarantee existence of minimizcrs for the cost wc are 
interested with, see [10] and Remark 6 below. An alternative model in which boundary 
conditions on directions are required with orientation is the one of Mumford. See 
Section 2.8.1. 
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Figure 2.1. A trajectory with two cusps. 



In this paper we are interested in the minimization of the fohowing cost, defined 
for smooth curves 7 in [6, c] : 



This cost is interesting for several reasons: 

• It depends on both length and curvature of 7. It is small for curves that are 
straight and short; 

• It is invariant by rototranslation (i.e. under the action of SE{2)) and by 
reparametrization of the curve, as should be any reasonable process of recon- 
struction of interrupted curves. 

• Minimizers for this cost do exist in the natural functional space in which this 
problem is formulated, without involving sophisticated functional spaces or 
curvatures that become measures. Indeed, in [10] we have proved: 
Proposition 5. For every {xb,yb),{^c,yc) e with {xb,yb) ^ (xcVc) and 
Vb, Vc £ K^\{0}, the cost (2.2) has a minimizer over the set 



Remark 6. In [34, 43, 44[, it has been proved that minimizers for the cost 
(2.2) are analytic functions for which 7 = at most for two isolated points. 
At these points the limit of ||7(i) ||iir^(t) is well defined. They are cusp points, 
i.e. points at which 7 becomes opposite. See Figure 2.1. 
Notice that at cusp points the limit direction (regardless of orientation) is 
well defined. In [10[ it is proved that if boundary conditions are required 
with orientation, then the cost (2.2) has no minimum over the set H^. 

However, the most interesting aspect from the modelling point of view is that this 
cost is a Riemannian length for lifts of planar curves over PTM^ (more precisely J [7] 
is a sub- Riemannian length, see below). In the spirit of the model by Petitot, Citti 
and Sarti. this is the most natural distance that one can define on PTM^. Indeed, this 
distance takes into account the fact that two configurations (xi, yi, 0i) and (x2, 2/2, ^2) 




(2.2) 
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are close each other if they are both close in the planar coordinates {x, y) and in the 
angle coordinate 6. 

Apparently, this cost is a good model to describe the energy necessary to excite the 
orientation columns that are not directly excited by the image (since they correspond 
to the corrupted part of the image) . Indeed it is a standard fact in sub-Ricmannian 
geometry (see Section 2.2) that the minimization of the cost J [7] is equivalent to the 
minimization of the energy-like cost 



The term ||7(i)|P models the energy necessary to activate horizontal connections, 
while the term ||7(t)||^i4'^(f) models the energy necessary to activate vertical connec- 
tions, see Figure 1.1. This is much more evident in the optimal control formulation of 
Section 2.3, where ||7(i)||^ is the control responsible for the "straight movements on 
M^" and ||7(t)pii'^(t) is the control corresponding for the "rotational movements on 
R^". Other models for these energies are of course possible, but our choice appears to 
be the most natural since it provides a well-posed variational problem. 

Finally, a key consequence of this choice of the cost is the following: we have a 
diffusion equation naturally associated with the variational problem. This diffusion 
equation can be used to reconstruct more complicated images than just curves. We 
use this diffusion as the key tool for the reconstruction algorithm. 

Remark 7. One could argue that there is no reason to give the same weight to 
the length term ||7|| and to the curvature term ||7(f)|piC^(t). However, if we define 
the cost 



with a fixed /3 7^ and if we consider an homothety [x, y) ^ {f3x, (3y) and the cor- 
responding transformation of a curve 7 = {x{t),y{f,)) to 7^ = {j3x{t), I3y{t)), then it 
is easy to prove that [7/3] = (3^J\ [7] = J [7]. Therefore the problem of mini- 
mizing is equivalent to the minimization of J with a suitable change of boundary 
conditions. 

Although the mathematical problem is equivalent by changing /3, this parameter 
will play a crucial role in the following, see Remark 21. 

Another interesting feature is the uniqueness of this sub-Riemannian distance. 
Beside the possibility of adding a weight (3 on the curvature term, that can be removed 
via an homothety, it is the unique sub-Riemannian distance for lift of planar curves 
on PTM^ that is invariant under the action of SE{2). See Proposition 14 below. 

2.2. Sub-Riemannian manifolds. In this section we recall some standard def- 
initions of sub-Riemannian geometry, that we use in the following. We start by re- 
calling the definition of sub-Riemannian manifold. 

Definition 8. A {n,m)- sub-Riemannian manifold is a triple (M, A,g), where 

• M is a connected smooth manifold of dimension n; 

• A is a smooth distribution of constant rank m <n satisfying the HormEinder 
condition, i.e. A is a smooth map that associates toq G M a m-dim subspace 
k{q) of TqM and \/ q <E M we have 





span{[Xi, [. . . [Xk-i,Xk] . . .]](g) | Xi e VecniM)} = T^M 
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where Vsch {M) denotes the set of horizontal smooth vector fields on M, 
i.e. 



VeCij(M) = {X e Vec(M) | X{q) G A{q) y q G M} . 

• gq is a Riemannian metric on A(g), that is smooth as function of q. 

A Lipschitz continuous curve q(-) : [0,T] — > M is said to be horizontal if q{t) (E 
A{q{t)) for almost every t € [0,T]. Given an horizontal curve q{-) : [0,T] — >■ M, the 
length of q{-) is 

«(?(•)) = ^smm,m) dt. (2.4) 

The distance induced by the sub- Riemannian structure on M is the function 

d{qo,qi) =mi{l{q{-)) \ q{0) = qo.q^T) = qi, q{-) horizontal}. 

The connectedness assumption for M and the Hormander condition guarantee 
the finiteness and the continuity of d(-, ■) with respect to the topology of M (Chow's 
Theorem, see for instance [6]). The function (i(-, •) is called the Carnot-Caratheodory 
distance and gives to M the structure of metric space (see [8, 23]). 

It is a standard fact that /(?(•)) is invariant under reparametrization of the curve 
q{-). On one side, if an admissible curve g(-) minimizes the so-called energy functional 

Eiqi-)) = f gq(t)im,m) dt. 
Jo 

with fixed T (and initial and final fixed points), then v = \/gg(t)(9(i)j <l{t)) is constant 
and q{-) is also a minimizer of /(•). On the other side, a minimizer q{-) of /(•) such 
that V is constant is a minimizer of E{-) with T = l{q{-))/v. 

A geodesic for the sub-Riemannian manifold is a curve q{-) : [0,T] — >■ M such 
that for each sufficiently small interval [^1,^2] C [0, T], then q{-)\n^ ^ minimizer 
of E{-). A geodesic for which gq{t){q{t),q{t)) is identically equal to one is said to be 
arclength parameterized. 

Locally, the pair (A,g) can be specified by assigning a set of m smooth vector 
fields spanning ▲, that are moreover orthonormal for g, i.e. 

A((z) = span {Xi(g),...,X^ ((?)}, g,(X,(q), X, (rj)) = (2.5) 

Such a set {Xi, . . . , X^} is called a local orthonormal frame for the sub-Riemannian 
structure. When (A,g) can be defined by m globally defined vector fields as in (2.5) 
we say that the sub-Riemannian manifold is trivializable. 

Given a (n, TO)-trivializable sub-Riemannian manifold, the problem of finding a 
curve minimizing the energy between two fixed points qo,qi £ M is naturally formu- 
lated as the following optimal control problem 

m 

g(t)=^«,(f)Xi(g(t)), (2.6) 

i=l 

Ui{.) e L°°([0,r],M), / ^u-(f) dt ~> min, 

q{0)=qo, q{T)=qr. (2.7) 
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It is a standard fact that this optimal control problem is equivalent to the minimum 
time problem with controls ui, . . . . Um satisfying ui{t)'^ + . . . + u„i(t)'^ < 1 in [0,T]. 
When the sub-Riemannian manifold is not trivializable, the equivalence with the 
optimal control problem (2.6)-(2.7) is just local. 

When the manifold is analytic and the orthonormal frame can be assigned by m 
analytic vector fields, we say that the sub-Riemannian manifold is analytic. In this 
paper we deal with an analytic sub-Riemannian manifold. 

A sub-Riemannian manifold is said to be of 3D contact type if n = 3, m = 2 and 
for every q £ M we have span{A((j'), [A, A](g)} = TqM. This is the case that we study 
in this paper. For details, see [5]. 

Remark 9. As a consequence of the invariancc by reparamctcrization of the cost 
(2.4), it is equivalent to state the minimization problem in the space of Lipschitz or 
absolutely continuous curves (i.e. for Ui(-) e L°°([0,T],R) or forUi(-) € L^{[0,T],M.).) 
See [10, Lemma 1]. 

2.2.1. Left-invariant sub-Riemannian manifolds. In this section we present 
a natural sub-Riemannian structure that can be defined on Lie groups. All along the 
paper, notations are adapted to group of matrices only. For general Lie groups, by gv 
with g in the Lie group G and v in the Lie algebra £, we mean {Lg)^,{v) where Lg is 
the left-translation on the group. 

Definition 10. Let G be a Lie group with Lie algebra £ and ^ C £ a subspace 
of £ satisfying the Lie bracket generating condition 

Lieq3:=span{[pi,[p2,...,[p„-i,Pn]]] | Pi € q3} = £. 

Endow <p with a positive definite quadratic form (.,.). Define a sub-Riemannian 
structure on G as follows: 

• the distribution is the left-invariant distribution 

A(5) :=ff*P; 

• the quadratic form g on A is given by 

Sg{vi,V2) := {g^^vi,g^'^V2). 

In this case we say that (G, A,g) is a left-invariant sub-Riemannian manifold. 

In the following we define a left-invariant sub-Riemannian manifold by choosing 
a set of m vectors {pi, . . . ,pm} which form an orthonormal basis for the subspace 
*P C £ with respect to the metric from Definition 10, i.e. ^ = span {pi, . . . , pm} and 
(pi,pj) = Sij. We thus have 

Hd) = gV = spa,n{gpi,...,gprn} 

and gg(s'pi,5pj) = Sij. Notice that every left-invariant sub-Riemannian manifold is 
trivializable. 

2.3. Lift of a curve on PTR^ and the sub-Riemannian problem. Consider 
a smooth planar curve 7 : [b, c] — > . This curve can be naturally lifted to a curve 
7 : [b,c] — >■ PTR^ in the following way. Let {x{t),y{t)) be the Euclidean coordinates 

of j{t). Then the coordinates of 7(t) are {x{t),y{t),9{t)), where 9{t) G M./{ttZ) is the 
direction of the vector {x{t), y{t)) measured with respect to the vector (1, 0). In other 
words, 

e{t) = arctan tt. (2.8) 



Of course we can extend by continuity the definition to points where 7(f) = but 
linif^f 9{t) is well defined. Wc assume 
[H] 9 : [b,c] R/(7rZ) is absolutely continuous. 

Notice that 6 = ||7||is'-y, hence hypothesis [H] is equivalent to require that ||7||if7 G 
L\[b,c],R). 

The requirement that a curve {x(t),y{t),9{t)) satisfies the constraint (2.8) under 
[H] can be slightly generalized by requiring that {x{t),y{t),9{t)) is an admissible 
trajectory of the control system on PTM?: 

± \ / cos{9) \ / 

y \ =ui{t){ sm{9) j+u2{t)lo | (2.9) 

with ui,U2 G L^{[b,c],R). Indeed each smooth trajectory 7 satisfying [H] is an 
admissible trajectory of (2.9). 

Since wi(t)2 = ||7(0IP, ^2(0' = = \\iit)fKy{t)^, we have 



Jb] = £ Vui{ty+U2{ty 



dt 



Hence, the problem of minimizing the cost (2.2) on the set of curves ^ is (slightly) 
generalized considering the optimal control problem (here q{-) = {x{-),y{-),9{-))) 

q = u,{t)X,{q)+U2it)X2{q), (2.10) 

/ cos(^^) \ / \ 

Xi{q) = I sm{9) \ , X^iq) = I J 1 ' (^.11) 

l{q{-))= I ^/ui{tf +U2{tfdt^mm, (2.12) 
Jb 

q{b) = {xb, Vb, 9b), q{c) = {xc, Vc, Oc), (2.13) 

{xb,yb) {xcVc), uuU2 e L\[b,c],R). (2.14) 

Remark 11. Notice that there are admissible trajectories g(-) = {x{-),y{-),9{-)) 

of the control system (2.10) for which the condition 9(t) = lim(_5.( arctan ^f^^ is not 

verified (consider for instance the trajectory x{t) = 0, y{t) = 0, 9{t) = t) or such that 
x{-) or y{-) fail to be smooth. However, it has been proved in [10] that minimizers of 
(2.10)-(2.14) are minimizers of (2.2) on the set ^ and they are smooth. 

Remark 12. (non-trivializability) A certain abuse of notation appears in 
formulas (2.9), (2.12), and (2.13), as in [45]. Indeed the vector field Xi is not well 
defined on PTR^. For instance, it takes two opposite values in 9 and 9 + n, that are 
identified. A correct definition of the sub-Riemannian structure requires two charts: 
• Chart A: 9 e]0 + kn,TT + knl k eZ, x,y eR. 



q = utit)X^{q) + U2{t)X2{q), X^ 
l{q{-)) = ^ut{tY+U2{tYdt, 
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• Chart B: 6* g] - 7r/2 + kn, tt/2 + knl k eZ, x,y eR. 



( cos(0) \ 

q = uf{t)X?{q)+U2{t)X2{q), = I sin(0) 1, 

l{q{-)) = £ ^uf{ty+U2{tydt, 

One can check that the two charts are compatible and that this sub-Rieniannian 
structure is non-triviahzable, while PTM? is parallelizable. 

Since the formal expression of Xi and Xi are the same, while they are defined 
on different domains, one can proceed with a single chart (however, one should bear 
in mind that ui changes sign when passing from the chart A to the chart B in ]R x 
Mx]7r/2, 7r[). In the following, since we study a "sum of squares" hypoelliptic diffusion 
on this sub-Riemannian structure, the problem disappears. 

This sub-Riemannian manifold is of 3D contact type: the distribution has dimen- 
sion 2 over a three-dimensional manifold and 

span{Xi(g), X2{q), [X^ X2M} = T^PTR^ . 



2.4. The sub-Riemannian problem on SE{2). It is convenient to lift the 
sub-Riemannian problem on PTM^ (2.10)-(2.14) on the group of rototranslation of 
the plane SE{2), in order to take advantage of the group structure. It is the group of 
matrices of the form 

cos{6) — sm{9) X 
SE{2) = { I sin(6i) cos(6i) y 
1 



9 G R/(27rZ) 
x,y gR 



In the following we often denote an element of SE{2) hy g = {x,y,9). 
A basis of the Lie algebra of SE{2) is {^1,^2,^3}, with 






Pi = u u u , P2 = i u u , P3 



We define a trivializable sub-Riemannian structure on SE{2) as presented in 
Section 2.2.1: consider the two left- invariant vector fields Xi{g) = gpi with i = 1,2 
and set 

kig) = span{Xi((;),X2((;)} gg{Xi{g),Xj{g)) = 6ij. 
In coordinates, the optimal control problem 

gGA{g), l{gi.))= f Jsg{t){9,9) dt ^ mm, (2.15) 

g{h) = {x\y\e% g{c) = {x'^,y-,e% (2.16) 
{x\y^)^{x'^,y% (2.17) 

has the form (2.10)-(2.14), but 6 G M/(27rZ). Notice that the vector field (cos(6l), sin(^), 0) 
is well defined on SE{2). 
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Figure 2.2. A case in which 6 £ M/(27rZ) is not the direction (with orientation) of'j. 



Remark 13. It is worth mentioning that the problem (2.15)-(2.17) (i.e. the prob- 
lem (2.10)-(2.14) with G K/ (27rZ)) cannot be interpreted as a problem of reconstruc- 
tion of planar curves where initial and final positions and initial and final direction of 
velocities (with orientation) are fixed. For instance, consider the curve starting from 
{x, y, 9) = (0, 0, 0) and corresponding to controls u\{t) = 7r/2 — t, U2{t) = 1. The cor- 
responding trajectory in the {x,y) plane is (— cos(t) + i(7r — 2f) sin(f) + I, tt sin^ (5) + 
tcos{t)—sm{t)). Notice that this trajectory has a cusp at time t = 7r/2. For t S [0, 7r/2[ 
we have that 6 is the angle with respect to (1,0) of the vector {x{t),y{t)), while for 
t e]7r/2,7r], it is not. See Figure 2.2. 

The control problem (2.10)-(2.14) defined on PTM? is left-equivariant under the 
action of SE{2). Indeed, topologically, PTR^ can be seen as the quotient of SE{2) by 
Z2 (in other words SE(2) is a double covering ofPTR2). In coordinates, {x,y,e) e 
PTB? corresponds to the two points {x, y, 9), {x, y,6 + n) € SE{2). Also, there is a 
natural transitive action of SE{2) on PTE? given by 

/ cos(6l) - sin(6') ( ^' \ ( cos{9)x - sm{9)y + x' \ 

sm{9) cos{9) y \ I y' \ ^ i sm{9)x + cosl9)y + y' j 

\ 1 / \ 9' J \ e' + e J 

eSE{2) ePTM2 ePTK2 

where ^' -|- ^ is intended modulo tt. The orthonormal frame for the sub-Riemannian 
structure on PTR^ given by Xi and X2 in formula (2.10) is indeed left-equivariant 
under the action of SE{2). 

In other words, given {x, y, 9) e PTK^ such that g e SE{2) satisfies {x, y, 9) = 
3(0,0,0), then 

Xi{x,y,e)=gpi, X2{x,y,9) = gp2. (2.18) 

The following proposition can be checked directly. 

Proposition 14. Let {PTM.^ , be a sub-Riemannian manifold and assume 

that it is left-equivariant under the natural action of SE{2). This means that if 
{Fi,F2} is an ortnonormal frame for the sub-Riemannian structure then 

F,{x,y,9)=gF,{0,0,0), F2{x,y,e) = gF2iO,0,0), (2.19) 

where g G SE{2) is such that {x, y, 9) = g{0, 0, 0) . Then, up to a change of coordinates 
and a rotation of the orthonormal frame, we have that 

( cos{9) \ / \ 

F,{x,y,9)={ sm{e) \ . F2{x,y,9) = { (2.20) 

V ; V1//3/ 
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for some (3 > 0. Notice that the problem of finding curves minimizing the length for 
the sub-Riemannian problem on PTR.'^ for which an orthonormal frame is given by 
(2.20), is equivalent to the optimal control problem (2.10), with the cost (2.4). 

2.5. The Sachkov synthesis. The solution of the minimization problem (2.10)- 
(2.14) on PTM^ can be obtained from that of the problem on SE{2) (2.15)-(2.17). 
The latter has been studied by Yuri Sachkov in a series of papers [34, 43, 44] (the first 
one in collaboration with I. Moiseev). 

The author computed the optimal synthesis for the problem. More precisely 
he computed the geodesies starting from the identity and for each geodesic the cut 
time, i.e. the time where it loses optimality. Thanks to the group structure, optimal 
geodesies starting from other points are just translation of these ones. In Figure 2.3 
the cut locus of the Sachkov synthesis is shown. 



Part of the cut locus due to topological reasons 




seen as a full torus with no boundary 

Figure 2.3. The cut locus of the Sachkov synthesis, i.e. the set of points where geodesies 
lose optimality for the sub-Riemannian problem on SE{2) (seen as the product of an open disc in 
R-^ times S^). Notice that the cut locus is adjacent to the starting point, as it always occurs in 
sub-Riemannian geometry. 

The complete optimal synthesis and the description of the cut locus for the prob- 
lem formulated on PTR^ has not been computed. However, as noticed by Sachkov, 
if we want to find the optimal trajectory joining {x,y,9) to {x,y,6) in PTM^, it 
is enough to find the shortest path among the four optimal trajectories joining the 
following points in SE{2): 

• {x,y,0) to {x,y,e) 

• (x, y,9 tt) to (x, y, 9) 

• (x, y, 9) to (x, y,9 + tt) 
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Figure 2.4. The problem of connecting level sets 



• (x, y,0 + tt) to (x, y,9 + tt) 
Moreover, Yuri Sachkov built a numerical algorithm for curve reconstruction on 
PTM^. In this paper, we will not go further on the subject of reconstruction of 
curves. For our purpose of image reconstruction, the sub-Riemannian structure only 
is important, since it allows to define intrinsically a nonisotropic diffusion process. 

2.6. The hypoelliptic heat kernel. When the image is not just a curve, one 
cannot use the algorithm described above in which curves are reconstructed by solv- 
ing a sub-Riemannian problem with fixed boundary conditions. Indeed, even if a 
corrupted image is thought as a set of interrupted curves (the level sets) , it is unclear 
how to connect the different components of the level set among them (see Figure 2.4). 

Moreover, if the corrupted part contains the neighborhood of a maximum or 
minimum, then certain level sets are completely missing and cannot be reconstructed. 

Remark 15. The difficulty of reconstructing a portion of an image contain- 
ing a maximum or a minimum is also the main drawbacks of methods based on 
sub-Riemannian minimal surfaces. These algorithms (see [15, 24, 45]) consider the 
boundary of the lift of the corrupted part as a closed curve 7 in the space SE{2) or 
PTR^. They then "fill the hole" with the surface that has boundary 7 and minimizes 
the surface area computed with respect to the sub-Riemannian metric. As clearly 
explained in [24], this method can fail for 3 main reasons: the minimal surface does 
not exist (depending on the regularity of 7), it can be non-unique, or even can exist 
but its projection on does not coincide with the corrupted part (either not covering 
the whole part or covering also a part of the non-corrupted image) . 

A second problem is that, even if the surface exists and it is computed, one has to 
choose how to diffuse the intensity of the image on the reconstructed surface. See ]24, 
Def 7.3] for the introduction of an "interpolation" function ft and a "disambiguation" 
function F. 

We then use the original image as the initial condition for the non-isotropic dif- 
fusion equation associated with the sub-Riemannian structure. This idea was first 
presented in [15]. 
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Roughly speaking, first we consider all possible admissible paths by replacing the 
controls in equation (2.10) by independent Wiener processes. Then, we consider the 
difi^usion equation which describes the density of probability of finding the system in 
the point {x, y, 9) at time t. 

More precisely, let {^i, .... X,„} be an orthonormal frame of a sub-Riemannian 
manifold and consider the stochastic differential equation 

m 

dqt = ^ X,{q)t dwl, 

i=l 

where the w"^ are independent Wiener processes. It is a standard result that, due to 
the Hormander condition, the stochastic process admits a probability density <j){q,t), 
that satisfies the Fokker-Planck-like diffusion PDE^ 

m 

dtcl>{q,t) = Y,X!cl>{q,t). (2.21) 

i=l 

For more details, see e.g. [31, 35]. 

Roughly speaking, the relation among the sub-Riemannian geodesies and the 
corresponding diffusion equation is the following: for small time the diffusion occurs 
mainly along optimal geodesies. 

For instance, a result due to Leandre [29, 30] states that if Pt{qi,q2) is the heat 
kernel associated to (2.21) then for f — ;> we have that —t\ogpt{qi, q2) — > d{qi, (?2)^/4, 
where d{-,-) is the Carnot-Caratheodory distance. For other results in this direction 
see [9, 27] and reference therein. 

In our case, this diffusion equation is 

dt(l){x,y,e,t)=AHcl){x,y,e,t) (2.22) 

where 

Ah = {Xif + {X2f = {cos{e)d^ + sm{e)dyf + dl 
Since at each point (a;, y, 0) we have 

span{Xi,X2, [Xi,X2]] = T(,^,,,,e)-P™', 

Hormander theorem [26] implies that the operator Ah is hypoelliptic. 

The diffusion described by the equation (2.22) is highly non isotropic. Indeed one 
can estimate the hypoelliptic heat kernel in terms of the sub-Riemannian distance 
(see for instance [4]), that is highly non isotropic as a consequence of the ball-box 
theorem (sec for instance [8[). 

Remark 16. Notice that the sub-elliptic diffusion equation corresponding to the 
sub-Riemannian structure (2.15)-(2.17) on SE{2), has the same form (2.22). The only 
difference is that on SE{2), e ]R/(27rZ). 

Remark 17. In [4] it has been proved that the Laplacian Ah is intrinsic on 
SE{2), meaning that it does not depend on the choice of the orthonormal frame for 
the sub-Riemannian structure. One can easily prove that this is the case also for Ah 
on PTM.^. 



^We emphasize here the fact that the PDE is not a stochastic one. 
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2.7. The hypoelliptic heat kernel on SE{2). The hypocUiptic heat kernel 
for the equation (2.22) on SE{2) was computed in [4, 17]. More precisely, thanks 
to the left-invariance of Xi and X2, the equation (2.22) admits a a right-convolution 
kernel pt{-), i-e. there exists pt such that 

e''^"M9) = ^o*Pt{9)= I Mh)Pt{h-'9Mh) (2.23) 

Jg 

is the solution for i > of (2.22) with initial condition 0(0, g) = (j>o{g) e V-{SE{2),M) 
with respect to the Haar measure fi. 

We have computed pt on SE{2) in [4]: 

/■+00 / +00 ■> 2 \2 

^'iCs) = / ^ E e""* < X)' ^^(5)cen(^, ^) > + 

+~ , \2 \ 

+ ^eM<se„(^,-),X^(ff)se„(^,-)> d\. (2.24) 

n=l / 

Here A indexes the unitary irreducible representations of the group and 

X\g) : L\S\C) L\S\C), 
X^{g)i}{a) ^ e^^(^ V (a + 0) 

is the representation of the group element g = (x, y, 0) on L'^{S^ , C). 

The functions se„ and the 27r-periodic Mathieu cosines and sines, and 

< 01,02 >:= Jgi (pi{a)(f>2{a.) da. The eigenvalues of the hypoelliptic Laplacian are 

:= — a„ (^x) ^-'^'-^ — ^" (^)' ^'^^'"^ '^^ ^^'^ ^" characteristic 

values for the Mathieu equation. For details about Mathieu functions see for instance 
[3, Chapter 20]. 

Since the operator dt — is hypoelliptic, then the kernel is a C°° function of 

{t,g) e R+ X G. Notice that ptig) = e^^^Suig)- 

The kernel (2.24) has been obtained by using the generalized Fourier transform. 
Once again, we refer to [4] for a detailed description of the generalized Fourier trans- 
form and the method to compute the kernel. 

2.8. The hypoelliptic heat kernel on PTR^. SE{2) is a double covering of 
PTR^. To a point {x, y, 9) £ PTM^ correspond the two points {x, y, 9) and {x, y, 9+^) 
in SE{2). From the next proposition it follows that we can interpret the hypoelliptic 
heat equation on PTR^ as the hypoelliptic heat equation on SE{2) with a symmetric 
initial condition. It permits also to compute the heat kernel on PTE? starting from 
the one on SE{2). 

Proposition 18. Let 0o S L^{SE{2). M) and assume that (l>o{x, y, 9) = (l>o{x, y, 9+ 
n) a.e. Then the solution at time t of the hypoelliptic heat equation (2.22) on SE(2), 
having 0o 0^ initial condition at time zero, satisfies 

<l){t,x,y,9)=^{t,x,y,9 + n). (2.25) 

Moreover if <po € L^(PTR^,R), then the solution at time t of the hypoelliptic heat 
equation on PTM^ (2.22) having 0o o,s initial condition at time zero is given by 

(l>{t,x,y,e) = [ (Po{x,y,9)Pt{x,y,e,x,y,9)dxdyd9 (2.26) 

JPTW^ 
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where 



Pt{x, y, e, X, y, e) := pt{{x, y, 9)-^ o {x, y, 9)) + pt{{x, y, 9)-^ o {x, y, 9 + 7r))(2.27) 



In the right hand side of equation (2.27), the group operations are intended in SE(2). 

Proof. Define the element 11 = (0, 0, tt) € SE{2) and observe the following prop- 
erties: 

• n is idempotent. 

• Property (2.25) reads as (f)Q{g\l) — <j)o{g)- 

• The kernel pt{g) satisfies ptiP-g) = pt{g^)- Indeed, call g = {x,y,9) and 
observe that, given a real function ip{a), we have 



Recalling the explicit expression of pt given in (2.24), we have ^((n^) = 
Pt {gll} . But Pt is real, hence the equality follows. 
We compute now ct){t,gir) in SE{2) with satisfying (2.25) and we prove that 
<t>{t-,g^) = 4>{t-,g)- Indeed, 



We now prove the expression (2.26) for 4){t, [g]) € i^(PTR^,M) for initial data 
4>o{[g\): where [g] is an element of PTK^, the class containing g and ^II in SE{2). 
Consider the function ^Joig) G L^{SE{2),'K.) defined by il^a{g) = (f'iig]): that clearly 
satisfies (2.25). Consider the unique solution tp{t,g) of the hypoelliptic equation 
(2.22), that is given by il>{t,g) = tpo *Pt{g)- Since ip{t,g) = ilj{t,g^), the function 
(f>{t, [g]) := ijj{t,g) is well defined. 

It remains to show that defined above is the solution of (2.22) on PTR"^. Indeed 
dt(p = dti) = AhiP- Since the vector fields defining both on SE{2) and PTR^ 
coincide, then the difi'ercntial operators Ah defined on SE{2) and PTR^ coincide, 
hence Ani^ = Ah<P- Thus (j) satisfies (2.22). Since (/)(0, [<?]) = 0o([5]), then </> is the 
(unique) solution. 

The explicit expression (2.26) is a direct consequence of the definition [g]) := 
tp{t,g) and of the explicit expression of tp given in (2.23). Indeed, 



The expression (2.26) is recovered by writing g = {x,y,9), h = {x,y,9) and recalling 
that gU = {x,y,9 + TT). D 



(n o g) V(a) = {{-X, -y, 9)) ^a) = 



= X^{{x,y,9 + TT))ij{a)=X^ (5on)V(a). 




</>(*, b])=^(t, 5) = / Mh)Pt{h-^g)dh= / Mh)Pt{h-^g)dh 
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2.8.1. Oriented vs. non-oriented approach. One of the key points of the 
algorithm presented in this paper is that directions are considered without orientation. 
As mentioned above, this choice is forced by well-posedness arguments when using the 
sub-Riemannian cost. Other approaches which consider directions with orientation 
are possible, but with a different cost. 

The most celebrated is the one due to Mumford [32], which in control language 
reads: 

/ cos{0) \ f ^ \ r r 

q= I sin(0) j +u{t) ( j . j {I + I3^u{tf) dt = j {1 + 0^ K^^^^)) dt ^ mm, 
q{b) = {xb,yb,0b), q{c) = {xc,yc,Oc), 

Here {x,y,9) € x S^. Notice that this variational problem has the same form as 
the one treated in this paper (i.e. (2.10)-(2.11) for the energy functional Jl^{ui{t)'^ + 
P^U2{t)^) dt), but forcing m to be one and taking 9 G instead of P^. In this way 
trajectories are parametrized by arclength for the Euclidean metric on M? and not for 
the sub-Riemannian one and, as a consequence, there are no cusps. Notice, however, 
that since the energy functional is not invariant by reparametrization, geodesies have 
a different expression. (These geodesies are "elastica" curves, see also [41, 42].) 

With a procedure similar to the one described in Section 2.6, one can naturally 
associate a diffusion equation to this model. Then, one gets a diffusion equation with 
drift, namely 

dt(t>ix, y, e, t) = (Xi + {X^f) 4>{x, y, 9, t) = {(cos{e)d^ + sm{9)dy + dj) <^{x, y, 9, t). 

The level set of the image can be oriented for instance on the left (or on the right) of the 
gradient of the initial condition. This choice is well defined when the initial condition 
is a Morse function. In practice, people consider both diffusions with positive and 
negative drift to have a more "symmetric" impainting. This approach was followed in 
[18, 19] for contour enhancement. 

Apparently, in the community, some researchers prefer Mumford's model, while 
others prefer the Petitot-Citti-Sarti model presented in this paper. Mumford's model 
has the advantage of not producing cusps (which are not observed in psychological 
experiments, see [38]), while the model presented in this paper has the advantage 
of treating horizontal and vertical connections at the same level and allows a more 
natural lift of the image. 

2.9. The mathematical algorithm. In this section we describe the main steps 
of the mathematical algorithm for image reconstruction. In the next section we give 
some guidelines for numerical implementation. 



STEP 1: Smoothing of Ic Assume that the grey level of a corrupted image is 
described by a function Ic : Vc '■= V'^ \^ [0, oo[. The set 17 represents the region 
where the image is corrupted. The subscript "c" means "corrupted". After making the 
convolution with a Gaussian of standard deviations Ox — o'y > 0^, we get a smooth 
function defined on M?, which is generically a Morse function: 

/c =Xc * G{cFx,<Jy). 

^Xc is considered to be zero outside X'c- Moreover we assume ax = to guarantee invariance 
by rototranslations of the algorithm. 
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We recall that a smooth function : — )• K is said to be Morse if it has only isolated 

critical points with nondogcncratc Hessian. Roughly speaking, a Morse function is a 
function whose level sets are locally like those of Figure 2.5. 





diffeomorphic to level sets 
of a linear function 



Maximum or minimum 



Saddle point 



Figure 2.5. Level sets of a Morse function 



STEP 2: The lift of /c : ^ M to a function : PTR^ _j. R 

This is made by associating to every point {x, y) of the direction G M/ (ttZ) 
of the level set of /c at the point {x,y). This direction is well defined only at points 
where V/c 7^ 0. At points where V/c = 0, we associate all possible directions (see 
Figure 2.6). More precisely, we define the lifted support Sf, associated with the 
function / as follows, 

Sf = {{x, y, ^) e R2 X s.t. Vfc{x, y) ■ (cos(^), sin(^)) = 0}, 

where the dot means the standard scalar product on M^. Let 11 : 5/ — )■ be the 

standard projection {x,y,6) <E Sf ^ {x,y) G M^. Notice that if Vfc{x,y) 7^ then 
Il^^{x,y) is a single point, while if Vfc{x,y) = then Il~^{x,y) = M/(7rZ). 

Let us study the set Sf, when fc is a Morse function. If (x, y) G is such that 
V/c (a;, 2/) 7^ and t/ is a small enough open neighborhood of {x, y), then the lift of Sf 
is an orientable manifold in U xP^. See Figure 2.7 A. If (a;, y) is an isolated maximum 
of fc, and J7 is a small enough open neighborhood of (x, y) having a level set of fc as 
boundary, then Sf is a Mobius strip in t/ x P^. See Figure 2.7 B. The same happens 
when (x, y) is an isolated minimum or saddle point of fc- Indeed we have: 

Proposition 19. If fc : ^ M. is a Morse function, then Sf is an embedded 
2-D submanifold ofR"^ X P^. 

Proof. Consider the surface S G SE{2) given by the equation 

gix, y, 9) := cos{e)d^Jc + sin{e)dyfc = 0. (2.28) 

If (x, y,9) € S then (.x, y,0 + tt) G 5 as well and 5 is a double covering of Sf. It is 
enough to show that S* is a surface. 

At points {x, y,9) S where V/c 7^ then dg is non-zero. Indeed dgg = would 
imply that the vector V/c 7^ is orthogonal to two non-zero vectors. At points where 
V/c = we have 
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Figure 2.6. Lift of an image with a maximum point. 




Figure 2.7. In Figure A, we draw the support of the lift of a linear function with nonvanishing 
gradient. Figure B presents the support of the lift of a function in a neighborhood of a maximum 
point. 

which cannot be zero since the Hessian Hf^ of fc is non-degenerate by the Morse 
assumption. □ 



STEP 3: Lift of fc to a distribution in x supported on Sf 

Consider the distribution on x P^: 

fc{x,y,0) = .fc{x,y)6{g) 

where 5{g) is the Dirac-deha distribution associated with g{x,y,9) :— cos{9)dxfc + 
sin{6)dyfc in the sense of [22, p. 222]. This distribution is supported on Sf and it is 
canonicahy defined by fc- Notice that this choice is not crucial and there are other 
possibihties. For example, in [45] the Dirac delta is replaced by a a power of the 
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cosine of the angle, centered on the angle 9. 

Remark 20. This step is formally necessary for the following reason. The surface 
Sf is 2D in a 3D manifold, hence the real function /c defined on it is vanishing a.e. as 
a function defined on PTK^. Thus the hypoelliptic evolution of /c (that is, the next 
STEP 4) produces a vanishing function. Multiplying /c by a Dirac delta is a natural 
way to obtain a nontrivial evolution. 



STEP 4: Hypoelliptic evolution 

Fix T > 0. Compute the solution at time T to the Cauchy problem, 

/ dt<k{x, y, 0, t) = {{cos{9)d, + sm{0)dy)' + P^dl)<k{x, y, 0, t) 

{ ^{x,y,0,O) = Mx,y,e). ^^■'^^> 

Remark 21. In the formula above, the Laplacian is given by + thus 

it depends on the fixed parameter /?. This means that we use evolution depending on 
the cost Jff rather than J. Tuning this parameter will provide better results of the 
reconstruction algorithm. 



STEP 5: Projecting down 

Compute the reconstructed image by choosing the maximum value on the fiber. 

fT{x,y) = max(j){x,y,9,T). 

Again other choices are possible for this projection. 

Remark 22. The algorithm depends on two parameters. The first is the time 
of the evolution T, the second is the relative weight /3 in formula (2.30). A variant of 
this algorithm consists of re-iterating the steps above for very short difi^usion times. 
This idea was already presented in [15] to build a minimal surface. 

Remark 23. One main feature of this algorithm is that it does not need the 
knowledge of the corrupted part. As a consequence the diffusion acts also in the 
noncorrupted region. The larger the diffusion time, the more modified image in the 
non-corrupted region. This is very visible comparing Figure 3.1 (small diffusion time) 
and Figure 3.2 (large diffusion time). This is one of the weak points of this completion 
process, that certainly takes place in the VI cortex as a low-level process. It is the 
counterpart of the global character of the method. 

However, due to the highly nonisotropic character of the diffusion, this effect is 
not too visible from a global point of view. 

Modifications of the algorithm which keep the original image unmodified are sug- 
gested in [15, 45], by admitting the diffusion in the corrupted part only. Also, in 
the standard PDE-based image processing inpainting algorithms, this problem disap- 
pears. Indeed, one solves a (stationary) elliptic-like problem on the corrupted part 
with Neumann boundary conditions, not an evolution equation. See e.g. [13]. 

3. Numerical implementation and results. First we present the main lines 
of the algorithm used in our simulations. 
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3.1. STEPS 2-3: Lift of an image. The formal definition of the lifted function 
is hard to realize numerically for two reasons: the discretization of the angle variable 
and the presence of a delta function. 

Both issues are solved changing the definition of the lifted function: 

fc{x,y,e) = fc{x,y)(j){Wfc{x,y),0), 

where (j){O,0) ^ l/(2e)V6l e ]R/(7rZ) and 61) = (\>^[exg{v) - 0) where (/)i(/3) is the 
TT-periodic function assuming the following values over the interval [0,7r]: 



l/(2£) if /3 £ [I -£,!+£] 

otherwise, 



for a fixed e > 0. 

Since the space is discretized, the non-zero values of /c are no longer defined over 
a set of null measure, hence the discretized hypoelliptic diffusion gives non vanishing 
function for all T > 0. Thus, it is not necessary to perform STEP 3. 

3.2. STEP 4: Hypoelliptic evolution. In this section we give the crucial 
ideas to compute efficiently the hypoelliptic evolution (2.30). Here (., .) is the scalar 
product in and is the rotation operator of angle Q. 

First of all, the main feature of the noncommutativc Fourier transform is to 
desintegrate the regular representation of SE{2). This was the main ingredient of the 
computation of the hypoelliptic heat kernel in [4]. Using the Fotirier transform again, 
th(^ hypoelliptic heat equation is transformed into a family of parabolic equations. 
These are more suitable for standard numerical methods. 

Roughly speaking, the non-commutative Fourier transform /(A) of the function 
f{x, y, 0) =: f{X, 6), for A G M^, is an operator meeting: 

[f{K)i^]{0)= [ [ /(X,a)^(a + 0)dae2-^<«-^^'^>dX=(7^)(i?_,A),(3.1) 

where *e is the convolution with respect to the angular variable and is the 2-D 
Fourier transform with respect to the spatial variables X = {x,y). 

Then it is natural to consider the Fourier transform with respect to X. Indeed, 
apply this transform u — >■ u to the initial value problem: 



that gives 



dtu = Ahu 
u{0,X,0) = f,{X,0), 



dtU = (i'^dju - ^■k'^{xcos{0) +ysui{0)Yu 
u{O,X,0) = %{X,e). 



(3.2) 



(3.3) 



Hence, for each point in the Fourier space, we have to solve an evolution equation 
with Mathieu right-hand term. 

It is easy to solve explicitly (3.3) over PTM^, i.e. with e M/tt. This simply 
divides the computation time by 4. 

This is the principle of the algorithm, which is massively parallelizable, since we 
can solve simultaneously the equation (3.3) at each point of the Fourier space. 
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3.3. Results of image reconstruction. In this section we provide results of 

image reconstruction using ttie algorithm presented above. For these examples, we 
have tuned the parameters /3, that is the relative weight, and T, the final time of 
evolution. 

Notice again that this algorithm processes the image globally and does not need 
the information about where the image is corrupted. The counterpart is that it 
modifies the non-corrupted part too. 

We present three results. 

• Figure 3.1 shows an image which is corrupted in a small piece of it. Then 
the diffusion can be applied for a rather small time avoiding an important 
diffusion effect in the noncorrupted part. 

• Figure 3.2 shows a strongly corrupted image. In this case a larger diffusion 
time is necessary to "inpaint" completely the corrupted part. The diffusion 
effect is clearly much more important. However in our opinion the result is 
surprisingly good. 

• The residual vertical and horizontal stripes on Figure 3.2 are not due to 

numerical discretization (they do not occur in Figure 3.1). They are the result 
of the diffusion of the original (white) grid. This is again a consequence of 
the fact that the diffusion process is global, as explained in Remark 23. In 
the spirit of global completion, this drawback is more or less unavoidable. 

• Due to pixelization of the image, one could think that corruption along the 
diagonal is the worst situation. Figure 3.3 show that this is not the case. 




Figure 3.1. Reconstruction of an image corrupted on a small portion. Here the diffusion is 
applied for a small time 

A Genericity of Morse properties of Gaussian convolution. In this ap- 
pendix, we prove that, generically, the convolution of a function over a bounded 
domain "D c with a Gaussian G is a Morse function. In particular, we first prove 
in Theorem 26 that the set of functions I G L^{T>, R) the convolution of which with a 
Gaussian is a Morse function is residuaF in (2?, R) . We then prove in Theorem 28 
that the set of functions I G (X>, R) such that I*G restricted to a compact K cM? 
is a Morse function is open and dense. 



^We recall that a subset of a topological space is residual when it is a countable intersection of 
open and dense sets. 
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Figure 3.2. Reconstruction of an image deeply corrupted. A larger time of diffusion is necessary 





Figure 3.3. Reconstruction of an image corrupted on the diagonal 



Definition 24. Let Z, Y be manifolds, F : Z Y a map and W C Y 
a submanifold. We say that F is transversal to W at z € Z, in symbols F rTT ^W, if, 
where y = F{z), either y or y €W and 

1. the inverse image {TzF)~^(TyW) splits and 

2. the image {TzF){TzZ) contains a closed complement to TyW in TyY . 

We say that F is transversal to W, in symbols F (TT W, if FJTi for every z G Z. 

We recall that a closed subspace of a Banach space E splits when there exists 
a closed subspace G such that E = F ® G. 

Remark 25. If E is Hilbert, then every closed subspace splits. See [2, Prop. 
2.1.15]. 

Theorem 26. Let V be a bounded domain of the plane . Fix ax, cry > 0. 
Consider the convolution map^ 

( L'^{V,R) C°°(IR2) 
\ I ^ I*G, 



is considered to be zero outside T). 
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where G is the Gaussian centred at (0, 0) 



1 -^-^ 
G{x,y):=- e ^ ^. 

LetX:= {I G L'^{V,U.) s.t. T{I) is a Morse function}. Then, X is residual in L'^{V, IS 
Proof. The proof relies on parametric transversality Theorems. The version we 
use is Abraham's formulation, see [1, Th. 19.1], recalled in the following. 

Theorem 27. Let A,X,Y be C" manifolds, p : A C""(X, F) a C" represen- 
tation, W C Y a submanifold, and evp : X x A ^ Y the evaluation map. Define 
Aw C A by Aw = {a E A \ paW\ W}. Assume that: 

1. X has a finite dimension n and W has a finite codimension q in Y, 

2. A and X are second countable, 

3. r > max {0, n — q}, 

4. evpmw. 

Then, Aw is residual in A. We apply Theorem 27 with A = L'^{'D,R), X = R^, 
r = 2. We choose F = x M x x and p the 2-jets of r(X), i.e., 

( A ^ C'-{X,Y) 

'''■{x ^ (Hi, n2, r(x), a.r(x), a,r(x), dl,T(i), dlyT{x), a^^rcx)) 

where 

111 : < . x and II2 : < . x 

{ {x,y) X { {x,y) y 

are the canonical projections. 

We fix 

W = {{x,y,a,pi,p2,qi,q2,q3) & Y s.t. {x,y) € M^, pi =p2 = 0, gigs - g| = O} . 

A function X G C^(M^) is a Morse function if and only if 

evp{x,y,l) = pi{x,y) = {x,y,T{I){x,y),dj:{I){x,y),dyT{I){x,y), 

tyn'^){x,y),dly\ 



dl,T{I){x, y), dlV{I){x, y), 9^ r(X)(a;, y)) 



does not belong to W for all {x,y) € R^. 

Remark that W is not a manifold. However, it is an algebraic set and hence it 
is a finite union of manifolds. In the following, we apply Theorem 27 as if W were a 
manifold, with the understanding that the Theorem is applied to each component. 

We now verify each of the conditions 1-4 in Theorem 27. Condition 1 holds with 
n = 2 and g > 3 for each component of W. Condition 2 holds, since A and X 
are separable metric spaces and hence second countable. Condition 3 holds for each 
component of W . 

Now we verify condition 4, that is the transversality condition eVp (Ti W . Fix 
x,y,X such that evp{X,{x,y)) G W. Condition 1 in Definition 24 holds because of 
Remark 25. We now verify condition 2 in Definition 24, where Z = x A. In the 
following, we prove that {T^x,y,i)^^p){T{x,y,x){^'^ x ^)) is the whole T'evp(x,j/,i)i^- The 
map T(^x,y,x)^^p has the following triangular form 



T'(x,y,I)eVp = 













1 


* 








Txevp{x,y,X) ^ 



(3.4) 
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We are left to prove that the tangent mapping Txevp{x, y,I) is surjective in M x x 
MP. for arbitrary fixed. After a suitable change of coordinate, we can assume 

that ax = (Ty = 1 and that (0, 0) G V. Let £ > such that V D Q := [—e, e] x [— £, e]. 
Define the function in L'^{V, M) 



61{x,y) 



Co + cix + C2y + cax" 



C4:Xy + C5y^ 



G{x-x,y-y) 



restricted to Q, and zero in r>\Q. The map p is linear in I, thus Txevp{x, y, X) [61] 
evp{x,y,6X). Consider the linear operator 



evp{x,y,dl) = 



( jQSX{x,y)G{x-x,y-y)dxdy \ 

/q SI{x, y)diG{x — x,y — y) dxdy 
Jq S^ix, y)d2G{x -x,y-y) dxdij 
Jq S^ix, y)df^G{x -x,y-y) dxdy 
Jq SX{x, y)dl2G{x -x,y-y) dxdy 
V lQSI{x,y)di2G{x-x,y-y)dxdy J 



as a function of the 6 variables (cq, . . . , C5), and consider the linear system evp(x, y, 61) = 
(a,pi,p2, 91, 92, 93), where (a,pi,p2, 9ij 92, 93) G 1^ is fixed. A direct computation 
shows that the determinant of the system is iQ^Q25a^a^ ^ ^' ^^^^ system always 
has a solution, i.e. Txevp{x,y,I) is surjective. 

By applying Theorem 27, we get Aw residual in A. We now prove that Aw = X. 
Since I G X implies evp{x, y,X) ^ W, then px fTT W, hence A D X. 

Now let us prove the inclusion A c X. Let X G A and fix {x, y) G R^. 

Nonintersection claim : px{x,y) ^ W. 

Proof of the claim. By contradiction, let 

w = evp{x,y,X) G W. 

Since px (Tl (x,y)W, then {T(x,y)Pi) ^{x,y)^) contains a closed complement to T^W 
inT^y. 

Observe that 

dim {T^x,y)Pi) {T(x,y)^^) < dim {T(x,y)^^) = 2 

and codivaT^W > 3, thus {T(x,y)Pi) {^{x,y)^^ cannot contain a closed complement 
to in T^F. A contradiction. 

By applying the claim for each (x, y) G M^, we get that px is a Morse function. □ 
Theorem 28. Let K he a compact subset ofM.^ with non-empty interior. Under 
the hypothesis of Theorem 26, the set \k := |x G L'^i'D, M) s.t. px]^ is a Morse function^ 

is open and dense in L^('D, R). 

Proof. Applying the openness of nonintersection Theorem [1, Th. 18.1] and 
using the nonintersection claim, we get that Xk is an open subset of L'^{V,R). Since 
Xk D X and X is dense, then the conclusion holds. □ 
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